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Abstract:  In  Monte  Carlo  and  molecular  dynamics  simulations,  million'  of  geometrical 
configurations  are  explored.  A  quick  calculation  of  the  energy  of  a  many-body  system  is 
crucial.  Of  all  the  computational  techniques  in  quantum  chemistry,  the  extended  Hiickel 
method  is  one  of  the  most  economic.  Although  simplistic,  it  often  gives  valuable 
qualitative  results  in  appropriately  selected  applications.  It  can  be  modified  to  give 
reasonable  geometries  of  molecules.  A  combination  of  molecular  dynamics  and  Monte 
Carlo  simulations  with  the  extended  Hiickel  method  should  be  useful  in  studying  dynamic 
phenomena  near  surfaces  and  in  the  solid  state.  We  have  taken  the  first  step  by  studying 
the  reconstruction  of  the  Si(lll)  (V3  x  V3)  R30°-Boron  surface  structure  using  a  modified 
extended  Hiickel  method  and  Monte  Carlo  simulation.  Our  simulated  annealing  suggests 
that  the  reconstruction  involves  a  three-atom  rotation  inside  a  silicon  matrix  instead  of  a 
direct  exchange  of  the  boron  atom  and  the  silicon  atom  under  it 


Introduction: 


Computational  algorithms  and  computer  machinery  have  advanced  to  such  a  stage 
that  simulations  of  dynamic  phenomena  of  practical  interest  are  feasible.  Tremendous 
insight  has  been  gained  from  these  numerical  experiments  Yet  in  certain  sense,  simulation 
science  is  still  in  its  infancy.  There  is  a  lot  of  'art*  in  this  rapidly  expanding  field.  This 
makes  things  challenging,  and  it  prevents  computer  simulation  from  degenerating  into  a 
louune.  There  is  no  substitute  for  excellent  chemical  and  physical  intuition.  One  of  the 
goals  of  this  work  is  to  share  the  joy  and  pain  of  computer  simulation  with  the  reader. 

The  quest  for  the  global  minimum  of  a  multi-dimensional  function  remains  one  of 
the  most  challenging  and  important  issues  in  science  and  technology.  The  major 
breakthroughs  in  optimization  occurred  in  the  fifties  and  sixties.  Yet  for  most  cases  of 
practical  interest,  the  standard  algorithms  are  only  capable  of  determining  the  local  minima. 
One  frequently  resorts  to  the  brute-force  approach  in  which  the  optimization  procedure  is 
executed  for  a  large  number  of  initial  configurations,  say  N.  As  N  -»  °o,  the  absolute 
minimum  is  located  with  certainty.  In  practice,  N  is  always  finite  and  the  global  minimum 
may  remain  illusive. 

Simulated  Annealing,  since  its  inception  in  1982,  has  emerged  as  a  promising  tool 
in  searching  for  the  global  minimum.1,2  It  has  solved  the  classical  travelling  salesman 
problem  (TSP),  involving  the  determination  of  the  shortest  cyclic  tour  for  a  salesman  to 
visit  a  large  number  of  cities  in  turn.  TSP  is  a  NP-complete  problem,  that  is,  the  CPU  time 
to  compute  the  exact  solution  scales  exponentially  with  the  number  of  cities.  Other  fruitful 
applications  include  the  design  of  large  scale  integrated  circuits,  image  processing,  code 
design,  neutral  network  theory,  refinement  of  X-ray  crystallographic  data,  etc.312 

Central  in  this  method  is  Nature's  own  mechanism  of  determining  its  most  stable 
state.  At  thermal  equilibrium,  the  probability  of  a  system  at  an  energy  E  above  the  ground 
state  is  characterized  by  the  Boltzmann  distribution  function: 
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p(E)  =  exp(-E/kT)/Z(D 


where  Z(T)  is  the  partition  function  —  a  normalization  factor  depending  on  the  temperature 
T,  and  k  is  the  Boltzmann  constant.  Even  at  low  temperatures,  there  is  a  finite  probability 
of  a  system  being  in  a  high  energy  state.  Correspondingly,  there  is  also  a  chance  to  escape 
from  a  local  minimum  in  search  of  the  global  one.  This  process  is  often  speeded  up  by 
annealing.  A  solid  is  heated  to  a  high  temperature  and  then  cooled  down  slowly.  High 
temperature  ensures  a  high  probability  of  escape;  hearing  provides  more  energy  for  the 
particles  of  a  solid  to  overcome  an  activation  barrier,  slow  cooling  allows  ample  time  for 
them  to  rearrange.  In  this  way,  strain  and  crystal  imperfection  are  minimized.  Strictly 
speaking.  Nature's  global  optimization  as  implemented  by  us  may  require  an  enormous 
amount  of  time.  Yet  all  real  and  numerical  experiments  have  to  be  completed  within  a 
reasonably  short  period.  Trapping  in  a  local  minimum  is  still  possible  and  the  process  of 
annealing  is  frequently  repeated  several  times  in  the  laboratory  and  so  it  has  to  be  in 
numerical  experiments  as  well. 

In  1953,  Metropolis  translated  these  idea  into  a  numerical  algorithm13  which 
eventually  forms  the  basis  of  the  simulated  annealing  procedure  proposed  by 
Kirkpatrick.14,15  Given  a  system  of  particles  with  an  energy  Ej,  a  new  configuration  is 
generated  by  displacing  each  particle  by  a  random  amount.  The  energy  of  the  new 
configuration  E2  is  computed.  If  the  new  configuration  is  more  stable  than  the  old  one, 
EpEj,  the  random  perturbation  is  continued  with  the  new  state.  If  AE=E2-E1>0,  the  new 
configuration  is  accepted  with  a  probability  given  by  the  Boltzmann  factor,  exp(-AE/kT). 
If  the  new  state  is  rejected,  the  random  perturbations  are  repeated  on  the  old  configuration. 
Following  this  Metropolis  criterion,  the  system  eventually  evolves  into  thermal  equilibrium. 
The  sequence  of  configurations  generated  in  this  process  defines  a  Markov  chain.  The 
length  of  this  chain,  L,  is  in  turn  characterized  by  the  number  of  successful  transitions. 
The  temperature  is  then  reduced  in  step,  typically  by  a  factor  of  0.75  to  0.90.  For  each 
step,  the  Markov  chain  should  be  long  enough  to  allow  the  system  to  approach  equilibrium. 
Provided  L  is  large  enough  and  the  decrement  in  temperature  for  each  step  is  small  (the 
essence  of  slow  cooling),  the  system  attains  its  ground  state  at  low  temperature.  The  final 
solution  obtained  by  simulated  annealing  should  not  depend  on  the  initial  configuration. 


2 


However,  this  method  is  computationally  much  more  intensive  than  ordinary  iterative 
algorithms.  In  actual  simulations,  it  has  been  found  that  after  cooling  to  a  reasonably  low 
temperature,  switching  to  a  conventional  minimization  algorithm  such  as  the  steepest 
descent  method  leads  to  a  great  reduction  in  computational  time.1,2  It  is  also  common 
practice  to  absorb  the  Boltzmann  constant  k  into  the  temperature  and  write  the  Boltzmann 
factor  as  exp(-AE/T).  Now  the  temperature  T  has  the  dimension  of  energy.  In  the  field  of 
combinatorial  optimization,  the  Boltzmann  factor  is  usually  written  as  exp(-E/C)  where  E  is 
the  cost  function  and  C  is  the  control  parameter. 

Surface  reactions  and  solid  state  reactions  are  of  utmost  technological  importance. 
They  are  of  direct  relevance  to  catalysis  and  corrosion.  Although  a  full  quantum 
mechanical  treatment  of  these  processes  is  usually  impossible,  simulations  of  biochemical 
phenomena  in  the  past  two  decades  have  proven  that  semi-empirical  and  even  empirical 
approaches  often  give  valuable  information. 16,1 7 

The  extended  Hiickel  method  has  been  providing  valuable  insight  into  bonding  in 
discrete  molecules  and  solid  state  compounds.18-21  It  is  also  useful  for  elucidating 
adsorbate-surface  interactions.22  This  approximate  molecular  orbital  calculation  is 
transparent  and  economic.  The  main  drawback  is  its  incapability  of  predicting  bond  length 
reliably.  A  new  version  of  this  method,  recently  proposed  by  Calzaferri,  seems  to  be  able 
to  predict  molecular  geometry  accurately.23  The  main  features  of  this  method  are  a 
distance-dependent  k  and  the  inclusion  of  two-body  repulsions.  (The  two-body 
electrostatic  interactions  were  originally  derived  by  applying  the  Hellmann-Feynman 
theorem.  Notice  that  Calzaferri's  version  is  actually  a  variant  of  the  atom  superposition  and 
electron  delocalization  (ASED-MO)  technique  first  put  forward  by  Anderson.24)  By 
combining  simulated  annealing  with  the  improved  extended  Hiickel  calculations,  we  hope 
to  provide  some  interesting  results.  The  parametrization  is  specified  in  the  Appendix. 

For  a  starter,  we  try  to  optimize  the  bond  length  of  H2  by  simulated  annealing. 
The  computed  bond  length  is  0.77  A.  The  experimental  value  is  0.74  A.  For  methane 
(CH4),  the  initial  C-H  bond  and  H-C-H  angle  are  arbitrarily  set  to  1.60  A  and  90°, 
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respectively.  The  optimized  bond  length  and  bond  angle  are  1. 19  A  and  109.5°,  compared 
to  the  experimental  values  of  1.10  A  and  109.5°,  respectively.  In  other  words,  the 
tetrahedral  geometry  has  been  regenerated  from  a  square-planar  configuration,  but  the 
calculated  C-H  distances  are  long. 

Stimulated  by  these  encouraging  results,  we  incorporated  the  above  two  features 
into  our  band  programs.  The  first  case  we  would  like  to  study  is  the  adsorption  of  a  boron 
atom  on  the  Si(l  1 1)  surface.  Due  to  its  promising  technological  properties,  this  system  has 
been  the  subject  of  intensive  investigation  recently.25'30  The  experimental  probes  applied 
include  scanning  tunneling  microscopy  (STM),  low  energy  electron  diffraction  (LEED)  and 
synchrotron  X-ray  diffraction,  while  the  density  functional  calculation  is  the  major 
theoretical  tool.  At  low  temperatures,  the  boron  adatom  sits  on  top  of  a  second-layer 
silicon  atom  and  is  bonded  to  three  other  silicon  atoms  of  the  top  layer,  as  indicated  in  the 
top  view  1.  (The  empty,  dashed  and  solid  circles  represent  the  top-layer  silicon  atoms. 
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second-layer  silicon  atoms  and  the  adsorbate,  respectively.)  The  coordination  number  of 
the  adsorbate  is  thus  four  and  this  adsorption  site  is  called  the  T4  site,  2.  Note  the  unusual, 
propellane-like  bonding  situation.  This  mode  of  adsorption  gives  rise  to  a  (V3  x  V3)  R30° 
LEED  pattern.  On  annealing  ( now  real,  not  in  the  computer! )  at  1000  °C,  the  boron  atom 
substitutes  for  a  second-layer  silicon  atom.  The  T4  site  is  in  turn  occupied  by  a  silicon 
atom.  The  new  silicon  adatom  is  directly  over  the  boron  atom,  3.  Now  the  boron  atom  is 


bonded  to  five  silicon  atoms  and  is  located  at  the  S5  site  (S  for  substitution).  This  is  the 
first  well-established  case  for  subsurface  substitutional  doping  of  a  semi-conductor.  It  is 
also  a  unique  situation  where  atoms  originally  adsorbed  on  top  of  a  surface  are  finally 
embedded  under  two  layers  of  the  host  material.  Introduction  of  point  defects  in  this 
surface  structure  leads  to  negative  differential  resistance  (NDR)  behavior,  that  is,  a  decrease 
in  electric  current  with  increasing  voltage.  NDR  is  the  crucial  property  that  allows  devices 
to  be  used  as  fast  switches,  oscillators,  and  frequency  locking  circuits.31,32  Notice  that  in 
the  famous  Si(l  11)  7x7  reconstruction,  the  silicon  adatoms  are  also  located  at  the  T4  sites. 
Finally,  the  whole  process  involves  breaking  and  reforming  of  many  Si-B  and  Si-Si  bonds. 
The  typical  bond  energy  for  a  covalent  single  bond  is  of  the  order  of  4  eV.  Hence,  the 
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establishment  of  the  mechanism  of  reconstruction  would  be  a  real  challenge.  Since 
simulated  annealing  is  basically  a  computer  simulation  of  the  physical  annealing  procedure, 
this  method  should  be  ideal  for  this  purpose. 
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Mode!  for  the  Si(lll)  Surface 


Silicon  adopts  the  diamond  structure.  Its  face-centered  cubic(fcc)  unit  cell  consists 
of  a  two-point  basis  at  0,0,0  and  1/4, 1/4, 1/4.  In  other  words,  the  silicon  structure  is  made 
up  of  two  interpenetrating  fee  lattices  displaced  relative  to  each  other  along  the  body- 
diagonal  of  the  cubic  unit  cell  by  one-quarter  of  the  length  of  the  diagonal.  Every  silicon 
atom  is  tetrahedrally  coordinated.  The  nearest  neighbor  distance  is  2.35  A.  The  hexagonal 
(111)  surface  is  the  close-packed  plane.  In  the  unreconstructed  surface,  there  is  one 
dangling  bond  for  every  (lxl)  surface  unit  cell.  There  are  in  turn  three  dangling  bonds  in 
the  Si(lll)  (V3  x  V 3)  R30°  unit  cell.  The  three  dangling  bonds  can  be  eliminated  by 
putting  a  boron  atom  at  the  T4  site.  Since  the  covalent  radius  of  boron  (0.90  A)  is  much 
smaller  than  that  of  silicon  (1.18  A),  adsorption  of  the  former  at  the  T4  site  would  induce 
substantial  amount  of  stress  at  the  silicon  surface.  In  fact,  assuming  a  completely 
undistoned  Si(l  11)  surface  and  putting  the  boron  atom  at  the  T4  site  such  that  it  is  2.4  A 
from  the  three  top-layer  silicon  atoms,  the  distance  between  the  boron  atom  and  the  fourth 
silicon  atom  would  be  about  1.6  A.  This  is  quite  short,  though  such  separations  are  found 
in  [1.1.1]  propellane.3336  Relief  of  surface  strain  has  been  proposed  to  be  the  main  reason 
for  the  instability  of  the  T4  site  relative  to  the  S5  site  25  27  Nevertheless,  the  boron  atom  is 
five-coordinated  at  the  S5  site.  This  may  introduce  even  more  strain  into  the  silicon  lattice. 
We  thus  have  some  reservation  about  the  above  proposal. 

In  a  compromise  between  speed  and  accuracy,  the  Si(  111)  surface  is  modelled  by  a 
five-layer  slab  of  atoms.  The  boron  atoms  are  put  on  one  side  of  the  two-dimensional  slab 
only.  In  all  calculations,  the  (V3  x  V3)  R30°  unit  cell  is  preserved. 

Although  the  exact  geometry  of  the  T4  site  is  unknown,  some  useful  information  is 
available.  When  the  boron  atom  occupies  the  S5  site,  only  its  five  nearest  neighbors 
labelled  in  3  show  significant  displacements  from  their  original  positions.  The  other  silicon 
atoms  are  essentially  stationary.  For  aluminum  and  gallium  adsorbed  at  the  T4  site  of 
Si(lll),  the  LEED  patterns  are  compatible  with  a  C3v  point  group  symmetry.37,38  The 
same  is  true  for  boron  at  the  S5  site.  We  therefore  optimized  the  positions  of  the  boron 
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atom  and  the  five  silicon  atoms  labelled  in  2  by  imposing  C3v  symmetry.  The  symmetry 
constraint  reduces  the  number  of  variables  from  eighteen  (3x6)  to  five.  The  optimized 
geometry  is  sketched  in  4,  which  also  serves  as  the  initial  configuration  for  our  computer 
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simulation.  The  C3  axis  is  perpendicular  to  the  surface,  passing  through  the  adsorbate. 
Notice  that  only  the  optimized  bond  lengths  are  given.  All  the  other  Si-Si  bonds  are  fixed 
at  2.35  A.  To  get  a  rough  idea  of  the  accuracy  of  our  optimization,  we  repeat  the 
minimization  procedure  for  the  case  of  aluminum  and  gallium.  The  computed  bond  lengths 
are  compared  with  the  experimental  findings  (in  parentheses)  in  5  and  6.  The  experimental 
error  is  about  0.1  to  0.2  A.  The  relatively  large  error  is  perhaps  enough  to  convince  the 
reader  that  determination  of  bond  lengths  at  the  surface  is  extremely  difficult  Despite  many 
years  of  intensive  research,  no  surface  structural  tools  can  give  results  as  reliable  as  X-ray 
diffraction  for  single  crystals.  The  R  (agreement)  factor  for  LEED  is  usually  higher  than 
15%.29>37’39  Yet,  dynamical  LEED  is  already  one  of  the  most  flexible  and  reliable  probes 
currently  available  for  complicated  surface  structure.  As  far  as  the  theory  goes,  a  complete 
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geometry  optimization  or  an  increase  in  the  thickness  of  the  silicon  slab  may  improve  the 
computed  bond  lengths.  (It  is  not  uncommon  to  use  a  slab  of  ten  layers  or  more  to  model 
the  Si(l  11)  surface.)  However,  the  focus  of  this  work  is  surface  dynamics,  not  the  statics. 
We  do  not  intend  to  invest  an  excessive  amount  of  computational  time  on  the  latter. 

The  silicon  slab,  together  with  the  boron  atom  at  the  T4  site,  is  then  annealed  in  the 
vacuum,  that  is,  there  is  no  boundary  in  the  z  direction  and  the  silicon  slab  is  periodically 
extended  in  the  x  and  y  direction.  The  coordinates  of  the  boron  atom  and  the  labelled 
silicon  atoms  in  2  are  varied  simultaneously  and  randomly.  The  positions  of  other  atoms 
are  fixed.  One  obvious  danger  of  this  procedure  is  that  some  of  the  atoms  may  leave  the 
slab  and  vaporize  into  the  vacuum.  Fortunately,  this  is  not  likely  if  the  Markov  chain  is 
long  enough.  Since  our  initial  geometry  corresponds  to  a  local  minimum,  all  standard 
optimization  algorithms  will  fail  to  effect  the  transformation  from  the  T4  site  to  the  S5  site. 
Simulated  annealing  is  the  only  option  available. 
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Implementation  of  Simulated  Annealing: 


To  specify  a  cooling  schedule  in  simulated  annealing,  four  parameters  have  to  be 

defined: 

a)  The  initial  temperature. 

b)  The  final  temperature. 

c)  The  length  of  the  Markov  chain. 

d)  The  cooling  factor. 

Its  is  perhaps  fair  to  state  that  a  universal,  efficient  cooling  schedule  for  global  minimization 
has  not  yet  appeared. 

The  choice  of  initial  temperature  T0  is  probably  the  most  important  part  of 
simulated  annealing.  If  TQ  is  too  high,  the  system  will  not  feel  the  attraction  to  the  global 
minimum.  Essentially,  it  just  wanders  randomly  on  the  potential  energy  surface.  An 
excessively  small  T0,  however,  would  lead  to  trapping  in  a  local  minimum.  To  ensure  a 
fast  approach  to  equilibrium,  Kirkpatrick  et  al.  suggested  that  T0  and  the  initial  stepsize 
should  be  chosen  such  that  no  less  than  80%  of  the  transitions  are  accepted.14  A  much 
lower  acceptance  probability  (25  to  50%),  however,  seems  to  be  sufficient  for  the  Lennard- 
Jones  or  the  Coulomb  potential.1,2  The  rule  of  thumb  for  Metropolis  sampling  is  that  the 
rejection  ratio  should  be  about  50%.  Yet  there  is  a  claim  that  an  acceptance  ratio  of  10%  is 
most  cost-efficient40 

In  view  of  these  controversies,  we  rely  on  our  chemical  intuition.  The 
experimental  annealing  temperature  is  around  1000  °C.  An  educated  guess  of  the  activation 
energy  of  the  reconstruction  would  be  one  to  two  electron-volts.  A  value  of  2.5  eV  for  T0 
should  be  sufficient.  From  the  molecular  simulations  mentioned  previously,  we  discovered 
that  an  acceptance  ratio  as  high  as  66%  seems  to  work  well.  We  therefore  adjust  our 
stepsize  to  attain  this  ratio.  The  initial  stepsize  turns  out  to  be  0.12  A. 
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Typically,  for  a  many-body  system  in  which  the  coordinates  of  X  atoms  are  varied 
simultaneously,  the  rule  of  thumb  is  to  perform  100X  trial  moves  for  every  temperature.1’2 
This  is  currently  beyond  our  computational  capability,  because  it  takes  about  36  hours  of 
CPU  time  on  a  Sun4  workstation  to  complete  a  calculation  for  600  configurations. 
Research  on  the  travelling  salesman  problem  (TSP)  indicates  that  instead  of  fixing  the 
number  of  trial  moves  for  each  temperature,  it  may  be  more  meaningful  to  focus  on  the 
number  of  successful  transitions,  that  is,  the  length  of  the  Markov  chain  L.  Indeed  there  is 
an  existing  algorithm  for  the  TSP,  setting  L  to  ION  where  N  is  the  number  of  cities.41  In 
our  case  18  atomic  coordinates  are  varied  simultaneously  and  randomly.  Thus,  L 
corresponds  to  180.  This  is  still  too  computationally  intensive  for  us,  so  we  reduce  L  to 
120.  Preliminary  simulations  suggested  that  there  is  a  tendency  for  the  silicon  atoms  to 
vaporize  into  the  vacuum  if  L  is  less  than  120. 

The  cooling  factor  is  chosen  to  be  0.875.  If  we  perform  the  simulated  annealing 
for  50  cycles,  the  final  temperature  would  be  3.6  meV.  This  temperature  should  be  low 
enough  for  all  practical  purposes.  (25  meV  is  roughly  equivalent  to  300  K.)  This, 
however,  also  corresponds  to  about  24  days  of  CPU  time  on  the  above-named  system. 
Fortunately,  after  the  second  cycle,  the  boron  atom  has  shifted  downward  over  a  distance 
of  1.6  A  and  bonds  to  a  third-layer  silicon  atom.  One  of  the  top-layer  silicon  atom  has 
moved  up  by  about  2.3  A  and  becomes  an  adatom.  The  adatom  in  turn  bonds  to  two  other 
silicon  atoms  below  to  form  a  bridge  trapping  the  boron  atom  within  the  silicon  matrix. 
The  silicon  atom  in  the  second  layer,  which  was  initially  2.3  A  under  the  boron  atom,  now 
shifts  above  the  latter  and  bonds  to  the  adatom.  The  major  activation  barrier  may  have  been 
overcome  and  there  may  be  no  need  to  maintain  the  system  at  such  a  high  temperature. 
Thus,  after  the  second  cycle,  we  anneal  the  whole  system  at  a  temperature  between  0.5  to 
0.2  eV  for  nine  more  cycles. 

At  the  end  of  the  eleventh  cycle,  we  obtain  all  the  essential  features  of  the  surface 
reconstruction.  The  boron  atom  has  shifted  downward  to  the  second  layer;  the  silicon 
lattice  starts  to  reform.  Hence  we  set  the  temperature  to  zero  and  continue  the  random 
search  for  the  minimum.  (Actually  we  have  performed  the  annealing  process  for  one  more 
cycle,  that  is,  about  10  more  hours  of  CPU  time.  However,  there  is  little  improvement  in 
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the  convergence  of  energy  in  the  twelfth  cycle.  That  is  another  reason  for  quenching  the 
system.)  At  the  end  of  our  simulation,  the  C3v  symmetry  is  regenerated  within  an  error  of 
0.016  A  in  the  atomic  coordinates.  The  boron  atom  has  settled  down  at  the  S5  site.  The 
S5  site  is  computed  to  be  more  stable  than  the  T4  site  by  about  2.5  eV. 

In  7,  we  compare  the  computed  bond  lengths  in  the  final  configuration  of  our 


computer  simulation  with  the  experimental  findings  obtained  by  LEED  (R  factor  = 
24.4%). 30  Since  there  is  a  slight  deviation  from  the  exact  C3v  symmetry,  the  computed 
bond  lengths  listed  in  this  figure  are  the  averages  over  the  three  bonds  supposed  to  be 
symmetry-equivalent.  In  Figure  1,  we  sketch  the  initial  configuration,  the  final 
configuration  at  the  second  cycle,  the  final  configuration  at  the  eleventh  cycle,  and  the 
configuration  at  the  end  of  our  simulation.  For  clarity,  only  the  top  and  side  views  of  the 
six  atoms  whose  positions  are  varied  are  shown.  Two  atoms  are  connected  by  a  line  if 
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their  separation  is  less  than  2.85  A. 


Figure  1  here 


In  summary,  our  computer  simulation  suggests  that  the  reconstruction  involves  a 
three-atom  rotation  inside  a  silicon  matrix  instead  of  a  simple  exchange  of  positions 
between  two  atoms.  Each  of  the  three  atoms  moves  over  a  distance  of  more  than  2.3  A. 
Hence,  on  the  extended  Hiickel  potential  energy  surface,  the  local  minimum  corresponding 
to  the  T4  site  is  really  far  from  the  deeper  minimum  of  the  S5  site.  Apart  from  our  limited 
computing  resources  which  confine  us  to  optimize  at  most  eighteen  atomic  coordinates,  and 
the  imposed  translational  symmetry  (an  intrinsic  restriction  of  our  band  programs),  our 
search  for  the  reaction  pathway  is  otherwise  unbiased.  There  is  no  additional  assumption 
concerning  the  mechanism  of  the  surface  reconstruction.  The  second  constraint  mentioned 
may  not  be  as  severe  as  it  seems,  because  experimentally  it  has  been  found  that  the  (^3x^3) 
unit  cell  is  preserved  when  boron  atom  is  located  at  the  T4  and  the  S5  site.  We  hope  that 
our  proposed  mechanism  will  be  tested  in  the  near  future  by  isotopic  labelling.  To  our 
knowledge,  this  is  the  first  study  of  surface  reconstruction  by  simulated  annealing.  There 
is,  however,  a  recent  study  of  a  simpler  surface  reconstruction  by  molecular  dynamics.42 

Density  functional  theory  (DFT)  is  ideal  for  extended  systems  and  large 
molecules.43,44  Although  computationally  intensive,  it  tends  to  give  excellent  results, 
especially  for  free-electron  systems.  Incidentally,  our  tight-binding  calculations  are  most 
likely  to  fail  in  such  systems.  Central  to  DFT  is  that  the  ground  state  energy  is  a  unique 
functional  of  the  electron  density.  Once  this  functional  can  be  determined  exactly,  there  is 
no  need  for  any  iteration  in  density  functional  calculations.  The  gain  in  CPU  time  would  be 
enormous.  Given  the  intensive  research  devoted  to  this  subject,  analytically  by  the 
coordinate  scaling45-46  and  numerically  by  quantum  Monte  Carlo  simulations,47  this  goal 
may  be  achieved,  at  least  partially,  in  the  near  future. 

We  believe  that  simulated  annealing,  together  with  molecular  dynamics  and  density 
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functional  theory,  will  become  the  main  tool  for  investigating  surface  dynamics,  such  as  the 
formation  and  reactivity  of  the  Si(l  1 1)  7x7  reconstruction;  catalytic  reactions  such  as  the 
Fischer-Tropsch  reaction;  or  corrosion  processes  such  as  the  interaction  of  oxygen  and 
water  with  metal  and  semiconductor  surfaces.  Simulated  annealing  can  also  be  applied  to 
study  dynamic  phenomena  in  the  solid  state,  such  as  the  diffusion  of  oxygen  in  high  Tc 
superconductors.  The  experimental  activation  energy  for  oxygen  diffusion  in  YBa^C^O-j 
is  about  0.4  eV.  Another  possibility  would  be  the  elucidation  of  the  mechanism  for  the 
formation  of  C^.  The  applications  of  computer  simulation  are  open-ended.  One  of  the 
foci  of  our  future  research  will  be  to  combine  density  functional  theory  and  computer 
simulation  to  study  dynamic  processes  in  condensed  phases,  and  the  mechanisms  for  the 
formation  of  clusters  such  as  C^,  C70,  Al6  and  Fe13. 
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Conclusion: 


We  have  investigated  the  mechanism  of  the  Si(l  1 1)-(V3  x  V3)  R30°-Boron  surface 
reconstruction  by  simulated  annealing,  using  a  modified  extended  Hiickel  method.  Our 
rough  computations  suggest  that  the  reconstruction  involves  a  three-atom  rotation  inside  a 
silicon  matrix,  as  shown  in  Figure  1.  Due  to  our  limited  computing  resources,  we  had  to 
impose  some  constraints  in  the  search  for  the  reaction  pathway.  The  dynamic  aspects  of 
surface  phenomena  are  very  difficult  to  study  either  theoretically  or  experimentally.  We 
would  be  more  than  rewarded  if  this  work  can  stimulate  further  research  in  that  direction. 
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Appendix 

All  calculations  are  performed  by  the  extended  Hiickel  method  within  the  tight- 
binding  formalism.48"51  The  distance-dependent  k 

k  =  1  +  (k  +  A2  -  A4K)exp(-6(r  -  d0)) 

where  k  =  0.80 

A  =  (HSi  -  HjjVfHjj  +  Hjj) 

8  =  0.13 

anH  the  two-body  repulsions  have  been  explained  in  a  previous  paper.23  The  extended 
Hiickel  parameters  are  listed  in  Table  1.  The  k  points  are  generated  according  to  the 
geometrical  method  of  Bflhm  and  Ramirez.52 


Table  1  here 
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Table  1.  Extended  Hiickel  Parameters. 


Atom 

Orbital 

Hj;  (eV) 

C 

H 

Is 

-13.6 

1.30 

C 

2s 

-21.4 

1.63 

2p 

-11.4 

1.63 

Si 

3s 

-17.3 

1.38 

3p 

-9.2 

1.38 

B 

2s 

-15.2 

1.30 

2p 

-8.5 

1.30 

AI 

3s 

-12.3 

1.17 

3p 

-6.5 

1.17 

Ga 

4s 

-14.5 

1.77 

4p 

-7.8 

1.55 

Acknowledgement: 


Y.  W.  would  like  to  dedicate  this  paper  to  Professor  Ruslan  Minyaev,  with 
gratitude  for  his  inspiration  and  encouragement.  Our  work  was  carried  out  with  the 
support  of  the  Office  of  Naval  Research.  Acknowledgment  is  also  made  to  Jane 
Jorgensen  and  Elisabeth  Fields  for  their  expert  drawings. 


18 


References: 


1 .  Wille,  L.  T.  Nature  1986, 324,  46. 

2.  Wille,  L.  T.  Chem.  Phys.  Lett.  1987, 133,  405. 

3.  Diswas,  R.;  Hamann,  D.  R.  Phys.  Rev.  B,  1986, 34,  895. 

4.  Carter,  H.  W.  Computer  1986, 19,  19. 

5.  Rothman,  D.  H.  Geophysics  1986, 51,  332. 

6.  Vecchi,  M.  P.;  Kirkpatrick  S.  IEEE  Trans,  on  Computer-Aided  Design  1983,  2, 
215. 

7.  Leong,  H.  W.;  Wong,  D.  F;  Liu,  C.  L.  Proc.  IEEE  Int.  Conference  on  Computer- 
Aided  Design,  Santa  Clara,  November  1985,  226. 

8.  Geman,  S.;  Geman,  D.  IEEE  Proc.  Pattern  Analysis  and  Machine  Intelligence, 
PAMI-6  1984,  721. 

9.  Wolberg,  G.;  Pavlidis,  T.  Pattern  Recognition  Lett.  1985, 3,  375. 

10.  Aarts,  E.  H.  L.;  Korst,  J.  H.  M.  Proc.  2nd  European  Simulation  Congress, 
Antwerp,  September  1986,  391. 

1 1.  Bemasconi,  J.  J.  Physique  1987,  48,  559. 

12.  Brunger,  A.  T.  J.  Mol.  Biol.  1988,  203,  803. 

13.  Metropolis,  N.;  Rosenbluth,  A.;  Rosenbluth,  M;  Teller,  A.;  Teller,  E.  J.  Chem. 
Phys.  1953,27,  1087. 

14.  Kirkpatrick,  S.;  Gelatt,  C.  D.  Jr.;  Vecchi,  M.  P.  IBM  Research  Report  RC  9355, 
1982. 

15.  Kirkpatrick,  S  ;  Gelatt,  C.  D.  Jr.;  Vecchi,  M.  P.  Science  1983, 220,  671. 

16.  McCammon,  J.  A.;  Harvey,  S.  C.  Dynamics  of  Proteins  and  Nuclei  acids, 
Cambridge,  1987. 

17.  Brooks,  C.  L.  IE;  Karplus,  M.;  Pettitt,  B.  M.  Protein:  a  Theoretical  Perspective  of 
Dynamics,  Structure  and  Thermodynamics,  Advances  in  Chemical  Physics,  vol.  71, 
Wiley,  1988. 

18.  Hoffmann,  R.  J.  Chem.  Phys.  1963,  39,  1397. 

19.  Hoffmann,  R.;  Lipscomb,  W.N.  J.  Chem.  Phys.  1962 ,37,  2872. 

20.  Ammeter,  J.  H.;  Biirgi,  H.-B.;  Thibeault,  J.  C.;  Hoffmann,  R.  J.  Am.  Chem.  Soc 


19 


1978, 100,  3686. 

21.  Whangbo,  M.-H.;  Hoffmann,  R.  J.  Am.  Chem.  Soc.  1978, 100,  6093. 

22.  Hoffmann,  R.  Solids  and  Surfaces:  A  Chemist's  View  of  Bonding  in  Extended 
Structures-,  VCH:  New  York,  1988  and  references  therein. 

23.  Calzaferri,  G.;  Forss,  L.;  Kamber,  I.  J.  Phys.  Chem.  1989,  93,  5366. 

24.  Anderson,  A.  B.  J.  Chem.  Phys.  1975,  62,  1187. 

25.  Headrick,  R.  L.;  Robinson,  I.  K.;  Vlieg,  E.;  Feldman,  L.  C.  Phys.  Rev.  Lett.  1989. 
63,  1253. 

26.  Bedrossian,  P;  Meade,  R.  D.;  Mortensen,  K.;  Chen,  D.  M.;  Golovchenko,  J.  A  ; 
Venderbilt,  D.  Phys.  Rev.  Lett.  1989,  63,  1257. 

27.  Lyo,  I.  W.;  Kaxiras,  E.;  Avouris,  Ph.  Phys.  Rev.  Lett.  1989,65,  1261. 

28.  Kaxiras,  E.;  Pandey,  K.  C.;  Himpsel,  F.  J.;  Tromp,  R.  M.  Phys.  Rev.  B  1990 .41. 
1262. 

29.  Huang,  H.;  Tong,  S.  Y.;  Quinn,  J.;  Jona,  F.  Phys.  Rev.  B  1990,  47,  3276. 

30.  Ma,  Y.;  Rowe,  J.  E.;  Chaban,  E.  E.;  Chen,  C.  T.;  Headrick,  R.  L.;  Meigs,  G.  M  : 
Modesti,  S.;  Sette,  F.  Phys.  Rev.  Lett.  1990,  65,  2173. 

31.  Lyo,  I.-W.;  Avouris,  Ph.  Science  1989,  245,  1369. 

32.  Bedrossian,  P.;  Chen,  D.  M.;  Mortensen,  K.;  Golovchenko,  J.  A.  Nature  1989,  542. 
258. 

33.  Wilberg,  K.  B.;  Walker,  F.  H.;  J.  Am.  Chem.  Soc.  1982,  104,  5239. 

34.  Honegger,  Ev.;  Huber,  H.;  Heilbronner,  E.;  Dailey,  W.  P.;  Wilberg,  K.  B.  J.  Am. 
Chem.  Soc.  1985, 107,  7173. 

35.  Wilberg,  K.  B.;  Dailey,  W.  P.;  Walker,  F.  H.;  Waddell,  S.  T.;  Crocker,  L.  S.;  Newton. 
M.  J.  Am.  Chem.  Soc.  1985, 107,  7247. 

36.  Wilberg,  K.  B.  Acc.  Chem.  Res.  1984,  17,  379. 

37.  a)  Northrup,  J.  E.  Phys.  Rev.  Lett.  1984,  53,  683 

b)  Huang,  H.;  Tong,  S.  Y.;  Yang,  W.  S.;  Shih,  H.  D.;  Jor.a,  F.  Phys.  Rev.  B  1990. 
42,  7485. 

38.  Kawazu,  A.;  Sakama,  H.  Phys.  Rev.  B  1988,  37,  2704. 

39.  Ogletree,  D.  F.;  Van  Hove,  M.  A.;  Somorjai,  G.  A.  Surface  Sci.  1986, 173,  351. 

40.  Wood,  W.  W.;  Jacobson,  J.  D.  Proceedings  of  the  Western  Joint  Computer  Conference 
(San  Francisco),  1959,  p.  261. 


20 


41.  Press,  W.  H.;  Flannery,  B.  P.;  Teukolsky,  S.  A.;  Vettrtling,  W.  T.  Numerical  Recipes, 
Cambridge,  1986,  p.331. 

42.  Ancilotto,  F.;  Andreoni,  Selloni,  A.;  Selloni,  A.;  Car,  R.;  Parrinello,  M.  Phys.  Rev 
Lett.  1990,65,  3148. 

43.  Trickey,  S.  B.,  Ed.;  Density  Functional  Theory  of  Many  Fermion  Systems, 

Advances  in  Quantum  Chemistry;  Academic,  1990,  vol.  21  and  references  therein. 

44  Dreizler,  R.  M.;  Gross,  E.  K.  U.  Density  Functional  Theory,  Springer- Verlag,  1990. 

45.  Levy,  M.;  Perdew,  J.  P.  Phys.  Rev.  A  1985,  32,  2010. 

46.  Levy,  M.;  Yang,  W.;  Parr,  R.  G.  J.  Chem.  Phys.  1985,  83,  2334. 

47.  Ceperley,  D.  M.;  Alder,  B.  J.  Phys.  Rev.  Lett.  1980, 45,  566. 

48.  Hoffmann,  R.  J.  Chem.  Phys.  1963,  39,  1397. 

49.  Hoffmann,  R.;  Lipscomb,  W.N.  /.  Chem.  Phys.  1962, 37,  2872. 

50.  Ammeter,  J.  H.;  Biirgi,  H.-B.;  Thibeault,  J.  C.;  Hoffmann,  R.  J.  Am.  Chem.  Soc. 
1978, 100,  3686. 

51.  Whangbo,  M.-H.;  Hoffmann,  R.  J.  Am.  Chem.  Soc.  1978, 100,  6093. 

52.  Ramirez  R.;  Btthm  M.C.  Int.  Quantum  Chem.  1986, 30,  391. 


21 


